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ABSTRACT 

We study statistical properties of stochastic variations in pulse arrival 
times, timing noise, in radio pulsars using a new analysis method ap- 
plied in the time domain. The method proceeds in two steps. First, 
we subtract low-frequency wander using a high-pass filter. Second, 
we calculate the discrete correlation function of the filtered data. As 
a complementary method for measuring correlations, we introduce 
a statistic that measures the dispersion of the data with respect to 
the data translated in time. The analysis methods presented here are 
robust and of general usefulness for studying arrival time variations 
over timescales approaching the average sampling interval. We apply 
these methods to timing data for 32 pulsars. In two radio pulsars, 
PSRs B1133+16 and B1933+16, we find that fluctuations in arrival 
times are correlated over timescales of 10 — 20 d with the distinct 
signature of a relaxation process. Though this relaxation response 
could be magnetospheric in origin, we argue that damping between 
the neutron star crust and interior liquid is a more likely explanation. 
Under this interpretation, our results provide the first evidence inde- 
pendent from pulsar spin glitches of differential rotation in neutron 
stars. PSR B0950+08, shows evidence for quasi-periodic oscillations 
that could be related to mode switching. 
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1 INTRODUCTION 



Radio pulsars are superb clocks. The regularity of the arrival times of the pulsed emission, upon barycen- 
tering, subtraction of proper motion, and gradual spin down, rivals atomic clocks in stability. Perhaps even 
more interesting are the imperfections of pulsars as clocks; variation in pulse arrival times, timing noise, is 
present at some level in all pulsars, and has remained unexplained since the discovery of pul sars in 1967 
Timin g no ise is a very complex pr ocess, and there appear to be many contributing factors; 



Hobbs et al 



(fcoioh and lCordes fc Downs! i|l985l ) lave shown that the complexity of timing noise cannot be explained by- 
high frequency random walks in the pulsar spin parameters (see, for example, Bovnton et al. 19721 : Cordesl 



nig_ 
LL98 



1980) . Possible contributors to tim ing noise include variable coupling between t he crust and t he liquid inte- 

199cl ), stochastic adjustments to the star's fi gure Jcordesll 19931), variations ir 



rior (jAlpar et al 



1986 



Jones 



the e xternal electromagnetic spin-down torque on the star ([Cheng 



1987a 



Urama et al. 



20101 ) . and time variability of the interstellar medium between the pulsar and Earth (jLiu et al 



2006; 



Lvne et al 



2011 



A neutron star is a dynamic system consisting of a rigid crust, a liquid interior, and an active magneto- 
sphere. Variations in the spin rate of the crust are the result of driving torques on the crust, filtered by the 
response of the system to those torques. Generally in noisy systems it is possible to determine properties of 
the response function without knowledge of the forcing function. For example, thermal fluctuations of the 
current in a circuit containing a resistor of resistance R in series with a capacitor of capacitance C can be 
used to measure the circuit's intrinsic decay time RC, independent of the processes that drive the current 
fluctuations. In this paper we develop techniques with which to measure properties of the response function 
of the neutron star system. 

Using timing noise to probe the neutron star system is an old idea. For a neutron star with a damped 
rotational mode associated with, for example, friction between the crust and a portion of the liquid interior, 
stochastic perturbations from rotational equilibrium by the noise process would never relax completely. In 
this case, decoupling by fluctuations at frequencies higher than the damping frequency t^ 1 would increase 
the spectral power for all frequencies above r^ 1 by a magnitude deter mined by the ratio of the moments 
of inertia of the crust and the liquid to which it is imperfectly coupled I Lamb et al. 197§h . Past studies of 



timing noise power spectra have not revealed structure beyond the power-law that is the hallmark of a nois 



process, that is, no deviations from rigid-body rotation have yet been detected (jBovnton fc Deeter 



Bovnton 



1981 



noisy 
19791 ; 



Bovnton et al.l 11984 ). The only evidence to date that neutron stars do not rotate rigidly 



comes from a very different p henomenon than timing noise: the glitches, sudden increases in spin rate (see, 



for example. iLvne et al.ll200ol ). whose occurrence and subse quent recovery have been attributed to variable 



coup ling between the crust and the liquid interior (see, e.g.. lAlpar et al 



2011 



198 



4; 



Link et al 



199 



Pizzochero 



Why neutron stars have not shown signatures of deviations from rigid-body rotation in their noise 
spectra has been an important open question in neutron star physics for over three decades. 

Since the early work on timing noise, the quantity and quality of data have increased. In the largest 



and most comprehensive study of timing noise to date, iHobbs et alJ (|2010l ) presented timing residuals of 
366 pulsars over nearly 40 years, some with nearly daily monitoring, mainly by the Lovell Telescope at 
Jodrell Bank. Here we revisit the issue of using timing noise to identify properties of neutron star response 
with these high-resolution data using methods applied in the time domain to 32 radio pulsars. In two 
radio pulsars, PSRs B1133+16 and B1933+16, we find that fluctuations in arrival times are correlated over 
timescales of 10 — 20 d with the distinct signature of a relaxation process. Though the relaxation response 
could be magnetospheric in origin, we argue that damping between the neutron star crust and interior liquid 
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is a more likely explanation. Under this interpretation, our results provide the first evidence independent 
from pulsar spin glitches of differential rotation in neutron stars. 

While this paper focuses on identification of signatures of the underlying physical processes responsible 
for pulsar timing fluctuations, the mathematical characterization of timing noise also plays a crucial role 
in efforts to detect gravitational waves in pulsar timing arrays such as the North American Nanohertz 
Observatory for Gravitational Waves NANOGrav ( |http: / /www.nanograv.org| ) , the Parkes Pulsar Timing 
Array PPTA ( http://www.atnf.csiro.au/research/pulsar/ppta/), and the European Pulsar Timing Ar- 
ray EPTA ( http://www.epta.eu.org). In particular, if the underlying power spectrum of timing noise is 
known, the signal-to-noise ratio for the gravitational wave background could be substantially improved 



(van Haasteren et al 



2009). 



In §2 we describe the data that we analyzed. In §3 we describe our analysis methods. In §4 we present 
results of the analysis for five radio pulsars. In §5 we discuss our results. The details of the analysis methods, 
examples of applications to simulated data sets, and tests of robustness are presented in the Appendix. 



2 DATA 



The Green Bank data were collected using a 25-m radio telescope at the National Radio Astronomy Obser- 
vatory, Green Bank, West Virginia. This telescope monitored the pulsars on a near daily basis from 1989 
through 1999. Observations followed a fixed daily schedule (adjusted a few times over the course of the 
program). Integration time depended on source strength, but was typically of order 40 minutes, and was 
divided into several subintegrations. For each of two linear polarizations, a filter bank produced total-power 
values for each of sixteen 1-MHz spectral channels across a band centered at 610 MHz. (The channels were 
not completely contiguous due to radio frequency interference considerations.). The two polarizations were 
balanced (using measured syste m noise), summed, and folded modulo the pulse period using the Princeton 
Mark 3 data acquisition system ijStinebring et al.lll992T ). This procedure produced a single pulse profile for 
each spectral channel in each subintegration. Off-line, the folded data profiles were cross-correlated with a 
high precision standard template to produce pulse times of arrival (TO As). The TOAs from all channels and 
all subintegrations of any given pulsar on any given day were averaged to produce a single effective daily 
TOA, w hich was used for th e present analysis. These TOAs were then analyzed with the tempo2 software 
package (jHobbs et al]l200a ) to obtain the timing residuals, using the JPL DE405 ephemeris. In each case, 
fits were made for position, proper motion, rotational frequency, and rotational frequency derivative. 

The Jodrell Bank data were collected using the 13-m Lovell Radio Telescope with nearly daily obser- 
vations from 1991 to 2008. Observations were made at 610 MHz, with a 4 MHz bandwith. 30-60 minutes 
of observations were made for each pulsar, depending on the flux density of each object, and then divided 
into 1-3 minute subintegrations. The subintegrations were averaged to produce a singl e profile, which wa s 
convolved in the time domain with the corresponding pulse templat e to produce TOAs (jHobbs et al.| [2004). 
These TOAs were then analyzed with the tempo2 software package ( Hobbs et al 1l2006h to obtain the timing 
residuals, using the JPL DE200 ephemeris. 
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Figure 1. Timing residuals in microseconds for four pulsars, showing typical long-period wander. The data are from 
the Jodrell Bank and Green Bank data archives (see Table 1). The residuals were obtained using the TEMP02 software 
package (Hobbs et al. 2006). 



3 ANALYSIS METHODS 

Timing noise appears as stochastic wander of the pulse arrival time, with respect to a deterministic spin- 
down model, over all timescales; examples are shown in Fig. [1] In most cases, the dominant effect is that 
of low-frequency variation over years. To study possible correlations over much shorter timescales, we work 
in the time domain and first subtract the slow wander of the timing residuals through application of a 
high-pass filter. We then perform two types of correlation analyses on these filtered (that is, "whitened") 
data. 



3.1 High-Pass Filtering 

We divide the time series into contiguous non-overlapping intervals of width W and calculate the average 
value of the residuals in each interval. Using unweighted least-squares fitting, we fit the average values 
with a cubic spline, and subtract the spline from the original time series to obtain the whitened residuals. 
More details and examples of this filtering method are given in the Appendix. For most data sets, a broad 
range in choices of the filter width W effectively removes long-period wander without introducing spurious 
correlations. 
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3.2 Discrete Correlation Function 

The data we analyze are unevenly sampled, and often contain large gaps. It is desirable to measure au- 
tocorrelations without having to bin the data, which would restrict our analysis to a timescale no shorter 
than the longest gap in the data and would entail a severe loss in time resolution. Unev enly- sampled data is 



readily handled with the discrete correlation function DCF(t) (Ede lson fc Kroliklll988f ). which we calculate 



as follows. For a set of arrival time residuals Sti measured at times ti, we construct the matrix 

a 2 

where St is the mean of the data set and a is its standard deviation (the choice of normalization is arbitrary). 
This matrix is calculated for all possible pairs (SU, 5tj), each of which is associated with the time difference 
Atij = tj — ti between the i-th and j-th measurements. (Here and in the following, St will denote residuals, 
and At will denote differences in measurement times). Suppose that there are M pairs that satisfy the 
condition 

T- Ar/2 sC AUj < t + At/2, (2) 

where At is the width of the sampling window. The discrete correlation function is the time average of eq. 
([T]) for the pairs that satisfy eq. ((2]), that is, 

DCF(r) = -^DCFy. (3) 

Points that do not fall in the sampling window do not contribute to DCF(r). The window size is chosen 
to maximize resolution without loss of statistical significance; in practice, At can be taken to be almost as 
small as the average sampling time. This procedure uses every data point, without any significant penalty 
in resolution due to occasional large gaps in the data. The values of r th at are used are binned in units of 



At. The standard uncertainty in the DCF(r) is (jEdelson fc Krolik 



1988 



^dcf(t) = jj-^ {J2 [DCF ij - DCF(r)] 2 } 1/2 . (4) 



3.2.1 Lagged Dispersion 

As another statistic with which to measure correlations, we introduce the lagged distribution function, 
LDF{<5t, t}, the distribution of fluctuation differences separated in time by a lag r. Correlations in the data 
can be studied with the lagged dispersion LD(r), which we define as the dispersion of the lagged distribution 
function. The LD has the following properties for a data set that is correlated over a timescale r c (and 
without a resonance). Because the time series will resemble itself to some extent upon time translation by 
times r < r c , the LDF{5t, r} will be relatively narrow, and hence LD(r) will be relatively small. For r > r c , 
the dispersion is larger because the data around time t + r are uncorrelated with the data around time t; 
in this case LDF{<5t,r} is broader than at low lag, and LD(r) asymptotes to some maximum value as r 
is increased. If the data series is simply noise, there will be no statistically-significant variations in LD(r) 
with r, since the data contain no timescale. A general increase of LD(r) up to a lag r ~ r c indicates that 
the series is correlated over a timescale ~ r c . Oscillations in the data appear also as oscillations in the LD. 
To handle uneven sampling, we construct LD(r) in a similar way as the DCF(t): 
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where Stij = 5ti — Stj. As before, the product is over the set of M elements that satisfies eq. ((2}, in which 
each element of the set is associated with the lag value r. The mean St in eq. © is defined as 

The LD is mathematically similar to the structure function used in studies of the turbulent interstellar 



plasma (e.g.. lRickettlll99dl : lYou et al.ll200'it l 



The standard uncertainty in LD(r) is 

^ld(t) = - Ttf - LD(t)] . (7) 

In the Appendix, we illustrate the usefulness of the DCF and the LD, in combination with high-pass 
filtering, to identify short-timescale correlations that are not readily identified with Fourier techniques. We 
find that high-pass filtering of the time series followed by calculation of the DCF or LD is a robust method 
for identifying an intrinsic relaxation timescale r c , provided the following conditions are satisfied: 

A^samp <C T c <C T wan j er , (8) 

where A£ samp is the mean sampling interval and r wan< j e r is the shortest timescale of the wander. If the first 
inequality is not satisfied, then the time resolution of the data is not sufficient to resolve the correlation 
timescale. If the second condition is not met, then the correlation cannot be disentangled from the wander. 
In practice, we must also require 

Tcorr < W < i wander i 

(9) 

to ensure that the filter removes the wander but not the correlation. In the extreme case of W ~ r sam p, the 
filtering method introduces spurious anti-correlations, indicating that the choice of W is too small. 



4 CORRELATION ANALYSES 

We have applied the methods of §3 to the 32 radio pulsars of Table [T] We apply a high-pass filter to 
the residuals, and then calculate the DCF and the LD. In most cases, we see no statistically significant 
correlations after the data are whitened. Here we present results for PSRs B1133+16, B1933+16, B0525+21, 
B1556-44, B0950+08, as particularly interesting examples that show correlated structure in their time series. 

PSRs B1133+16 and B1933+16. These two pulsars have been monitored almost daily with the 
12.8 m telescope at Jodrell Bank, usually at 610 MHz. Timing residuals for PSR B1133+16 over about 4000 
days, and PSR B1933+16 over about 6000 days, are shown in Fig. [2] We begin with these pulsars because 
they show particularly interesting correlations in their times series. 

In Figs. [2] and [3] (3rd panel) we show the DCFs for B1133+16 and B1933+16, after whitening with 
filter widths W of 400 d and 120 d, respectively. Uncertainties were obtained from eq. (J4j . (All uncertainties 
in this paper are 1-a). The DCF for B1133+16 shows highly significant, non-periodic correlations over a 
timescale r c ~ 10 d, with the distinct structure of a relaxation process (see Appendix). B1933+16 shows 
similar correlations, over a timescale r c ~ 20 d. From the calculated uncertainties, shown in the figures, the 
confidence level of the detected correlations is high. We estimate the correlation significance to be 
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Table 1. The pulsars to which we have applied the analysis techniques of this paper. JB - Jodrell Bank timing data, 
from the 13m dish with an observing frequency of 610 MHz, with a 4 MHz bandwidth, using 30-60 min observing 
sessions. GB - Green Bank timing data, from telescope 85-3 (25m diameter), with an observing frequency of 610 
MHz and 16 MHz bandwith, with typical observation sessions of ~ 40 min. The pulsar parameters given are from 
ATNF Pulsar Catalog (http://www.atnf.csiro.au/people/pulsar/psrcat/). 

where the product is over all values of n < r c , the inferred correlation timescale. The quantity 5* should 
give an approximate measure of the significance of the result in units of the (typical) statistical uncertainty 
of each data point. For B1133+16, the correlation over r c = 10 is seen over the first four bins of width 
At — 2 d, giving S = 24.5. For B1933+16, the correlation over r c = 20 is seen over the first nines bins of 
width At — 2 d, giving S = 31.2. These values of S indicate very high statistical significance. As a check, 
we obtain a lower limit on the significance with a bootstrapping method. We create an ensemble of data 
sets in which the time-ordering of the whitened residuals is randomly shuffled. The DCF is calculated for 
the shuffled data, with A = 2r d for both pulsars. We calculate the significance S a huffled for r ^ r c (the 
first four bins for PSR B1133+16, and the first nine bins for PSR B1933+16). We then ask: what fraction 
of shuffled sets give Scuffled > S? Under the null hypothesis of uncorrelated data, Scuffled is comparable to 
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Figure 2. Timing fluctuations in PSR B1133+16, in microseconds. Second panel - timing fluctuations after sub- 
traction of the long-term wander evident in the top panel, using a filter width W = 400 days. The uncertainties are 
comparable in magnitude to the fluctuations, and are not shown for clarity. The gaps in the data sets are on account 
of equipment upgrades; our analysis techniques enable us to use data on both sides of the gaps. Third panel - DCF 
of the whitened residuals. Fourth panel - LD of the whitened residuals. Bottom panel - power spectrum of the 1700 
day span of unwhitened residuals (see text). 



S, since S represents one particular realization. In 10 4 shufflings for each data set, the null result was never 
realized. We therefore estimate the statistical significance of the detected correlations to exceed 1 — 10 -4 , 
consistent with the values of S obtained above. This estimation method has the advantage that nothing is 
assumed about the underlying statistics of the data, rather, the data themselves are used to evaluate the 
likelihood of the null hypothesis. 

In Figs. and [3] we show LD(t) for PSRs B1133+16 and B1933+16. The value of LD for uncorrelated 
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Figure 3. Same as Fig. [2] for PSR B1933+16. The filter width is W = 120 d, for which the correlations over 
timescales of ~ 20 d are most clearly seen. The power spectrum in the bottom panel is for the 1400 day span of 
unwhitened residuals (see text). 



data is calculated by averaging many shufflings of the data; this value is used to normalize LD(r) shown 
in the figures (see Appendix A3). For B1133+16, we find that LD(r) is significantly lower than what we 
expect from the null hypothesis out to r ~ 10 d, in support of the DCF result. B1933+16 is very similar, 
but indicates a relaxation timescale r c ~ 30 d, consistent with the DCF. 

The power spectra of the unwhitened data do not show the correlations that our time-domain analysis 
has revealed. Gaps in observations of 4-5 days are frequent for B1133+16, and several gaps of 8-16 days are 
present. To obtain a power spectrum without interpolating data, we have chosen a span of data approxi- 
mately 1700 days in length, with gaps no longer than four days. The power spectrum of the timing residuals 
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are shown in Fig. [2] (bottom panel) , using uniformly weighted bins of width four days. The spectrum is 
generally white, with some excess power at low frequencies. No high frequency features or spectral breaks 
are evident. For PSR B1933+16, we select a segment of data spanning 1400 days, also placed in bins of 
width four days, to calculate the power spectrum shown in Fig. [3] (bottom panel). The red spectrum corre- 
sponds to the low-frequency wander visible in the timing residuals (Fig. [3} . The Fourier transform spreads 
the power over all frequencies, obscuring correlations that are more readily identified in the time domain. 

The two complementary and independent statistics we have used here, DCF(r) and LD(r), measure 
different properties of the data. The DCF measures the extent to which the data set compares to itself upon 
translation in time, while the LD measures the extent to which the distribution of fluctuation differences 
widens as the data decorrelate for larger lags. These two different statistics give consistent results for PSRs 
B1133+16 and B1933+16. Several tests of the robustness of these results are given in the Appendix. In 
particular, the whitening process cannot introduce spurious correlations at timescales shorter than the 
window width. The filter width of W = 400 d used for the analysis of PSR B1133+16 is a factor of ~ 40 
larger than the derived correlation time. The filter width of W — 120 d used for the analysis of PSR 
B1933+16 is a factor of ~ 6 larger than the correlation time. 

PSR B0525+21. Timing residuals are shown in Fig. [4] (top panel). The wander in the timing noise 
occurs over timescales > 300 days. Upon whitening with a filter width W = 100 days, the DCF (bottom 
panel) shows variability over timescales of ~ 10 d, but no evidence of relaxation. Relaxation response, if 
present, could be overwhelmed by the variability. 

PSR B1556-44. Timing residuals are shown in Fig. [5] (top panel). DCFs for this pulsar are shown 
in Fig. [5] for W = 100 (middle panel) and 50 days (bottom panel). For W = 100 days, the DCF shows 
variability over ~ 10 d. With W = 50 days, more wander is removed, but variability over timescales > 10 d 
is still evident (bottom panel). As for PSR B0525+21, any possible signature of relaxation response might 
be overwhelmed by the variability. 

PSR B0950+08. Timing residuals are shown in Fig. [6] Upon whitening the data, the DCF of this 
pulsar shows both correlations and anti-correlations (Fig. [6] bottom panel), associated with quasi-periodic 
variability over a timescale of < 100 d. Th e pulse profile of PSR B095 0+08 changes modes over timescales 
comparable to those identified in the DCF jshabanova fc ShitovllioO^ ; the variability in the DCF could be 
directly related to mode changing. 



5 DISCUSSION 

We find strong evidence for correlation of timing residuals over timescales less than ~ 10 d in PSR Bl 133+16, 
and ~ 20 d in PSR B1933+16, indicative of relaxation response in the neutron star system of crust, 
core, and magnetosphere. While it is possible that the correlations are created by variable torques in the 
magnetosphere, we consider this possibility to be unlikely. The basic timescale in the magnetosphere is 
the light travel time across the light cylinder, which is comparable to the spin period. For the observed 
correlations to be magnetospheric in origin, the magnetosphere mus t have a re l axatio n time of weeks. 
Such a long relaxation time is difficult to explain on physical grounds. iLvne et al. I fcoioh have shown that 
many stars exhibit changes in spin-down rate in association with pulse-shape changes, clear evidence that 
the magnetosphere is dynamic, but these changes always occur much more quickly than ~ 10 d. PSR 
B1133+16 and PSR B1933+16 show no evidence for such modes changes, though PSR 1133+16 does show 
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Figure 4. Top panel - Timing residuals for PSR B0525+21, in microseconds. Middle panel - residuals after whitening 
with W=100 days. Bottom panel - The DCF of the whitened data, showing variability over ~ 10 d. 



nulling. Mode changes would show both correlations and anti-correlations, unlike B1133+16 and B1933+16, 
which show only correlations. 

The correlations cannot be due to time variability of the interstellar medium. The delay in the arrival 
time of a pulse of frequency i^mhz due to dispersion in the interstellar medium is (|Lvne fc Graham-Smith 
20051 ') 



At IS M = 4.15 x 10 6 MHz 2 pc _1 cm 3 ms x i^L x DM > ( n ) 

where DM is the dispersion measure (cm -3 pc). The dispersion measure is calculated during multi- 
wavelength observations for 131133+16 and B1933+16, and used to subtract the delay as part of the 
timing solution. The time derivative of the dispersion measure has been measured from long-term moni- 
toring: d(DM)/dt ~ 8 x 10" 4 cm" 3 pcyr" 1 for B1133+16, and ~ 2.3 x 10" 3 cm" 3 pc yr" 1 for B1933+16 
( Hobbs et al.ll2004l ). This variability of the electron density produces small fluctuations in the arrival time 
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Figure 5. Top panel - Timing residuals for PSR B1556-44. Middle panel - The DCF after whitening with a filter 
width W = 100 d. Variability over > 10 d is evident. Bottom panel - The DCF with W = 50 days. More wander has 
been removed, but variability remains. 



delay, on the order of 10 /isyr . To estimate the maximum effect of variations of the interstellar medium, 
we superimposed a sine wave of amplitude 10 /is with period 10 d on our whitened data. No detectable 
signal is introduced into either the DCF or the LD. 

An alternative explanation is that the correlations found for PSRs B1133+16 and B1933+16 are due 
to damping between the neutron star crust and interior liquid as the system is excited, presumably often, 
away from a state of rotational equilibrium that is never reached. In this picture, our results provide the 
first evidence independent from pulsar spin glitches of differential rotation in neutron stars. In support of 
this interpretation, we note that the measu red correlation ti mes are comparable to the post-glitch recovery 
timescales measured in many pulsars (e.g., Lvne et al. 200oh . 



The methods introduced here are also generally useful for identifying quasi-periodic processes that occur 
over timescales of days. PSR B0950+08 is a particularly clear example, showing altern ating correlations and 
anti-correlations in the DCF. This star exhibits mode changes over a similar timescale (|Shabanova fc Shitov 
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Figure 6. Top panel - Timing residuals for PSR B0950+08, in microseconds. Middle panel - The residuals after 
whitening with a filter width W = 100 d. Bottom panel - T he DCF of the whitened r esiduals, showing variability 
over timescales comparable to observed pulse mode changes (Shah anova fc ShitovlEx)! ). 



2004), and so the DCF appears to be showing magnetospheric effects. Though not as strong as PSR 



B0950+08, PSRs B0525+21 and B1556-44 also show evidence for variability. 

For most of the pulsars in Table [TJ we find no evidence of relaxation behavior. In order to resolve 
correlation timescales of ~ 10 d, nearly daily sampling is required. Most of the pulsars that we have 
analyzed do not have the high sampling rates of PSRs B1133+16 and B1933+16. These pulsars with lower 
sampling rates might have short-timescale correlations that are not detectable yet with the methods used 
here. For data which contain relaxation and periodicity over similar timescales, it is not possible to remove 
the wander without removing any non-periodic correlations that may exist. 

The analy sis methods applied a re suitable for seeking time correlated structure in any noisy time 
series (see, e.g. 



Fukumura et al 



2010 for a recent application of autocorrelation methods to the analysis of 



astronomical data). For pulsar timing noise in particular, these techniques show promise for distinguishing 
effects that might be magnetospheric in origin, such as mode switching, from those effects that are directly 
related to rotational response of the stellar interior. Mode switchin g, for example, i s now known to be 



related to sudden changes in magnetospheric torque in many pulsars (jLvne et al 



20101 ). This effect can be 



seen as quasi-periodic behavior of the DCF, as we have shown with PSR B0950+08. By contrast relaxation 
response, as seen in PSRs B1133+16 and B1933+16, is more naturally attributed to internal dynamical 
degrees of freedom. An in-depth study of relaxation response in the context of a physical model is in 
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progress. The techniques developed in this paper may also be used to detect precession and oscillation 
modes in pulsar timing data. 
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APPENDIX A: EXAMPLES AND TESTS 

In this Appendix we describe the analysis techniques used in this paper in more detail and present many 
examples and tests of robustness. Synthetic data sets are generated to simulate the properties of real radio 
pulsar timing data. We also discuss a simple physical model for random spindown torques that captures 
the essential features of the detected correlations. 



Al Whitening with a High-Pass Filter 

To illustrate the effects of the high-pass filter that we use, we construct simulated timing residuals 5t(t) 
consisting of periodic functions, quasi-periodic wander and gaussian noise, 

St(t) = sin M + <t>i) J + p (t) + N(t), (Al) 

where Ai, uji, and (j>i are the amplitude, frequency, and phase, respectively, P(t) is an arbitrary polynomial 
function, and N(t) is gaussian noise. Ai and 4>i are randomly generated values. The sine waves are weighted 
by the frequency uii to produce a red power spectrum, similar to that of typical pulsar data. An example 
of a simulated time series is shown in Fig. IA1I (top panel), using a sum of 10 sine waves. For this example, 
we choose frequencies uii — LJi/i, where i is an integer, and a weighting of wr . The maximum frequency is 
lji = 0.32 (period=19.6), and and the minimum frequency is ujio — 6.25 x 10~ 4 (period=10 4 ). 

We now whiten the synthetic data as follows. We divide the time series into contiguous non-overlapping 
intervals of width W and calculate the average value of the residuals in each interval. Using unweighted 
least-squares fitting, we fit the data in each interval with a cubic spline and subtract it from the original 
time series to obtain the whitened residuals. The effect of whitening for different values of the filter width 
W can be seen by calculating the autocorrelation function of the whitened residuals. The autocorrelation 
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Figure Al. Top panel - Simulated residuals (arbitrary units) consisting of 10 sine waves, an arbitrary 5th order 
polynomial, and Gaussian noise, and the corresponding ACF. Second panel - Whitened data with a filter width 
W = 300. Third panel- Further whitening with filter width W = 150. Bottom panel- Further whitening with filter 
width W = 75. 



function measures the similarity of a time series to itself upon translation in time by a given lag t, 

ACF(r) = Y, ~ M ^f t+r ~ M) , ( A2) 

i 

where Xt is the measured signal at time t, fi is the mean of Xt, and o~ is the square root of the variance of 
Xt- In this example, for W = 300 the polynomial has been removed well, but significant periodicity remains 
in the time series (Fig. IA1I 2nd panel). For W = 150, some of the periodicity has been removed, but high 
frequency components are still apparent (Fig. I All 3rd panel). Using a filtering parameter W = 75 removes 
almost all of the wander in the time series, leaving a correlation function consistent with Gaussian noise 
(Fig. I All 4th panel). Recall that the highest-frequency sine wave in the series has a period of 19.6, less than 
W, but the high-frequency components have been suppressed in proportion to the period. 

Anti-correlations arise if W is too small, comparable to the sampling interval. In this case, the cubic 
spline begins to fit the noise, introducing spurious power at high frequencies. This results in a whitened 
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Figure A2. Autocorrelation function of the simulated time series in Fig. IA1I (top panel), after whitening with 
W = 20, showing anti-correlation at low lag, the signature of overwhitcning. The anti-correlation is strongest for lags 
~ W/2. 



time series which is anti-correlated for lags r < W (Fig. IA2|) , as fitting the noise produces residuals which 
are more likely to be of opposite sign. To see this, consider the limiting case in which each bin of width W 
contains only two points. The average value of each bin is subtracted from the time series, resulting in data 
that are anta-correlated for a lag r = 1. In general, the anti-correlation due to overwhitening is strongest 
for a lag r ~ W/2. This distinct signature allows us to easily determine if the selected filter width W is too 
small. As a general rule, we need to keep W at least several times larger than the average sampling interval. 



A2 Comparison of the Time-Domain Analysis with Fourier Methods 

To illustrate the advantages of time domain versus frequency domain analysis for the type of time series 
that we analyze, we consider simulated data produced by a simple physical model for a neutron star with 
an internal degree of freedom. The interaction between the fluid component and crust component of the 
star is described by 

I c Q c (t) = N(t) - — (Q c (t) - Q f (t)) , (A3) 

/,%(*) = -0M*) (A4) 

where I c and If are the moments of inertia for the crust and fluid, respectively, Q c (t) and are the 

rotation rates of the crust and fluid, r c is the coupling timescale between the crust and fluid, and N(t) 
is the torque on the crust, from internal or external sources. To construct a simple model of neutron star 
spin behavior, we consider crust and fluid components initially in co-rotation, with a ratio of moments of 
inertial I c /If = 1 for illustration. We perturb the crust with a series of S- functions, 

N(t) = ^2Aid(t-ti), (A5) 

i 

where Ai are randomly generated amplitudes. This model for the torque represents a series of instantaneous 
transfers of angular momentum to the crust, regardless of the source or sources. In between impulses, eqs. 
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Figure A3. Top left - Simulated rotation frequency residuals with data correlated over a timescale r = 1. Top 
right - The power spectrum of the frequency residuals, indicating a "knee" at lot c ~ 1. Bottom left - The simulated 
residuals with Gaussian noise. Bottom right - The power spectrum of simulated residuals. After the addition of noise 
to the time series, the coupling timescale is no longer evident. 



31) and (IA4I) have the solution 



t/T 



(A6) 



where fii and ^2 are constants determined by Ai of the last 5-function, I = If + I c is the total moment 
of inertia, and r = (7/// c )r c . Using this model, we construct a time series of frequency residuals shown in 
Fig. IA3I consisting of 2000 points with an impulse applied every unit time, and a sampling rate of 10 points 
per unit time. 

The rapidly-perturbed system has relaxation response over a timescale r c , as can be seen from the 
response function in frequency. Fourier transforming eqs. ()A3[) and (IA4|) . and eliminating the Fourier trans- 
form of Qf(t), gives 



17 c (oj) = 



{0JT c f + (Ic/If)' 2 



\N{lo) 



(A7) 



(cor c y + {I/I S Y 

We define the response function as the frequency dependent term in brackets, shown in Fig. IA4I At frequen- 
cies higher than ut c — 1, only the crust responds to the torque. At lower frequencies, the entire moment 
of inertia of the star / is perturbed by the torque, resulting in a lower response. There is a "knee" in the 
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Figure A4. Frequency dependent response function of the crust in the two-component model, normalized to unity. 
The knee at ujt c represents the transition from coupled to decoupled response. 

spectrum at oir c ~ 1. Such a knee is evident in the power spectrum of the simulated phase residuals (Fig. 
IA3I - top right panel). 

In Fig. IA3I (bottom panels), we have added gaussian noise to the time series to simulate noisy data. 
In this case, the knee in the spectrum corresponding to the coupling timescale r c becomes buried beneath 
the noise. Time domain analysis is better suited for this noise-dominated time series. We calculate the 
autocorrelation function of the time series shown in Fig. I A3I (bottom left). The coupling time is readily seen 
in the autocorrelation function shown in Fig. IA5I 

A3 DCF Tests 

To illustrate that the DCF can distinguish between correlations due to wander and those due to a relaxation 
process, we add a series of 10 "sawtooth" functions to the synthetic data of Fig. IA3I (lower panel). Each 
sawtooth consists of a discontinuous jump of random magnitude, followed by a linear decay over 10 time 
units. The jumps are randomly spaced over the time series, with magnitudes drawn from a Gaussian with 
a standard deviation of unity about zero; see Fig. IA6I (top panel). The resulting data are shown in Fig. 
IA6I (bottom panel); the effects of the sawtooth are too small to be seen in these simulated data. We have 
also added gaps to simulate real data. We show data with even sampling (apart from the gaps); we have 
confirmed that uneven sampling has a negligible effect on the results. 

DCFs for this time series upon whitening are shown in Fig. IA7l for three filter widths W. For W = 75, 
only positive correlations over ~ 10 time units are evident (top panel), as the low- frequency wander has 
been almost completely subtracted. This DCF signature is easily differentiated from the variability arising 
from the impulsive noise, which produces both correlations and anti-correlations. At lower values of W, the 
fitting function begins to remove correlations in the data, and anti-correlations begin to appear (middle and 
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Figure A5. The autocorrelation function for the simulated time series shown in Fig. [A3] (bottom left panel), showing 
the intrinsic relaxation time of the system that is not visible in the power spectrum. 



bottom panels); these data have been overwhitened. For a broad range in W, correlations remain evident, 
even if the data have been overwhitened. 

The LD for the synthetic data of Fig. IA6I (bottom panel) is shown in Fig. IA8I We decorrelated the data 
through random shuffling of the points. For 100 such shufflings, we calculated the average LD(r) to establish 
the base LD(r) one would expect under the null hypothesis that the data are completely uncorrelated. These 
values are shown as the horizontal lines in each figure; LD(r) is normalized in terms of this average value. 
This method has the advantage that nothing is assumed about the underlying statistics of the data, rather, 
the data themselves are used to evaluate LD(r) for the null hypothesis. Uncertainties are calculated using 
eq. [Jj The LD shows that the data are correlated over a timescale of 10 time units, consistent with results 
of the DCF for this simulated time series (Fig. IA7j ). 

We conclude that high-pass filtering of the time series followed by calculation of the DCF is a robust 
method for identifying an intrinsic relaxation timescale t c , provided the following conditions are met: 

^^sainp ^ 7"c ^ ^wander ; (^^) 

where At sam p is the mean sampling interval and r wan der is the shortest timescale of the wander. If the first 
inequality is not satisfied, then the time resolution of the data is not sufficient to resolve the correlation 
timescale. If the second condition is not met, then the correlation cannot be disentangled from the wander. 
In practice, we must also require 

< w < ^wander ) 

(A9) 

to ensure that the filter removes the wander but not the correlation. In the extreme case of W ~ r samp , the 
filtering method introduces spurious anti-correlations, indicating that W is too small. 
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Figure A6. Top - Schematic of the sawtooth function. Each impulse decays linearly over a timescale of 10 units. 
Impulses are randomly spaced in time, with random amplitudes of either sign drawn from a Gaussian distribution 
with a standard deviation of unity. Bottom - Simulated residuals consisting of periodic functions, a polynomial, 
Gaussian noise, and a sawtooth function. Several gaps are added to simulate the sampling we have for the real data. 
The effects of the sawtooth function are too small to be seen against the wander and added noise. 

A4 Robustness Tests 

To verify that the correlations found using the DCF are not introduced by the whitening method, we 
compare our results to those found using other whitening methods, using data from PSR B1133+16 for 
illustration. The wander in this pulsar can be adequately removed with a fifth-order polynomial. Applying 
the DCF analysis to the data whitened using this technique, the results are very similar to those we find 
using the technique described in Section 2.1 fFig. I A9l middle panel). We also apply the DCF to unwhitened 
data for PSR B1133+16, and then whiten the correlation function, rather than the time series, to remove the 
signature of low- frequency wander from the DCF fFig. IA9l bottom panel). These three whitening techniques 
give very similar results. 

For B1133+16, the large difference between the relaxation timescale (~ 10 d) and the wander timescale 
(~ 1000 d) allows a wide range of choices for the high-pass filter width W that give similar results. The 
wander is not well removed for W > 600 days, and the fitting function begins to subtract the relaxation 
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Figure A7. DCFs for whitened simulated residuals using W = 75, 50, 25, after adding a "sawtooth" function to 
the simulated residuals shown in Fig. IA6| The relaxation timescale of 10 units is easily identified for W = 75, with 
low-frequency wander nearly completely subtracted by the high-pass filter. For W = 50, the fitting function begins 
to subtract the correlations. At W = 25, the time series is overwhitened, resulting in reduction of the relaxation 
signature and anticorrelations from lags in the range 5 < r < 25. 



correlations for W < 100 days. We show DCFs for PSR B1133+16 for several values of W in Fig. [KM Cleft 
column). 

The characteristic timescale of the wander is much shorter for 131933+16, about 300 d. We plot several 
DCFS for B1933+16 for different values of W (Fis. lATOl right column). To successfully remove the wander, 
W must be a factor of several smaller than the shortest characteristic timescales of the wander. We also 
ensure that W is a factor of 2-3 greater than the relaxation timescale that appears in the DCF as W is 
reduced. In the case of B1933+16, these constraints leave little room to vary W. For W ~ 80 — 120 days, 
the DCFs are nearly identical. For Figs. [5]and[31 we use a cutoff frequency of / ~ 400 d _1 for B1133+16, 
and / ~ 120 d" 1 for B1933+16. 
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Figure A8. LD of the simulated time series shown in Fig. lA6l after whitening with W = 75. The relaxation process 
occurring over 10 time units is evident, confirming the results of the DCF shown in Fig. IA7I 
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Figure A9. DCF for PSR B1133+16, using three different whitening methods. In the top panel, the whitening 
method described in Section 2.1 is used with W = 400 d. (For other values of W, see Fig. lAlOt . For the middle 
panel, timing residuals were whitened using a fifth-order polynomial. For the bottom panel, the DCF function itself 
was whitened using W = 400 d (see text). 
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Figure A10. DCFs for PSR B1133+16 (left column) and PSR B1933+16 (right column) for several values of W. 
For B1133+16, the top panel uses W = 600 days, the middle panel uses W = 200 days, and the bottom panel uses 
W = 80 days. In the range W ~ 100 — 500 days, the DCFs are nearly indistinguishable. For B1933+16, the top panel 
uses W = 140 days, the middle panel uses W = 100 days, and the bottom panel uses W = 60 days. For W ~ 80 — 120 
days, the DCFs are similar and show correlations for r < 30 d. 
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